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On  a Method  of  Numerical  Differentiation 


Mira  Bozzini  and  Milvia  Rossini 


Abstract.  In  this  paper  we  present  a method  for  the  numerical  differen- 
tiation of  two-dimensional  functions  when  scattered  data  are  given.  The 
method  is  based  on  a regularization  of  the  given  sample. 


§1.  Introduction 

In  this  paper  we  present  a method  for  the  numerical  differentiation  of  two- 
dimensional  functions  when  scattered  data  are  given.  The  problem  of  numer- 
ical differentiation  is  very  important  when  dealing  with  function  approxima- 
tion. In  fact,  a satisfactory  recovery  of  a function  given  in  a sampled  form 
needs  some  knowledge  of  the  derivatives. 

It  is  also  well  known  that  this  problem  is  ill-conditioned.  Consider  for 
instance,  one-dimensional  equispaced  data  with  hi  = 10'  and  h2  = 2P.  On  a 
computer,  because  of  the  base  change,  numerical  differentiation  gives  consid- 
erable errors  in  the  first  case,  and  a more  accurate  solution  in  the  second  case. 
Moreover,  it  is  strongly  influenced  by  the  data  position. 

The  literature  on  scattered  data,  includes  the  papers  [7,8,10],  their  im- 
provements [3,5],  and  some  experiments  on  their  use  [9].  These  papers  provide 
the  gradient  approximation  at  the  sampled  points  (see  [7,10])  or  the  approxi- 
mation of  the  gradient  function  (see  [8]).  This  is  done  by  triangulation  or  local 
and  global  moving  least  square  interpolation.  In  some  of  them,  asymptotic 
bounds  for  the  error  are  also  supplied. 

Generally,  these  methods  provide  a satisfactory  approximation  inside  the 
domain  in  which  the  data  are  given,  but  they  may  give  large  errors  at  the 
boundary  (see  Figs.  1-4  below). 

In  this  paper  we  present  a method  based  on  a regularization  of  the  sample 
which  gives  an  error  with  an  uniform  behaviour  on  the  whole  domain  in  which 
the  data  are  assigned.  The  regularization  is  done  by  constructing  a new  set 
on  a regular  grid.  Namely,  taking  into  account  the  previous  observations,  we 
have  considered  dyadic  grids.  Obviously  the  construction  of  this  new  set  will 
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generate  errors  that  can  be  thought  of  as  causal  errors  depending  on  the  sam- 
ple. Then  (see  [4]),  to  smooth  the  data,  we  perform  a wavelet  decomposition 
of  the  signal.  Using  the  smoothed  data,  we  approximate  the  gradient  at  the 
grid  points  by  the  classical  finite  centered  difference  formulas,  and  finally  we 
construct  a smooth  function  approximating  the  unknown  gradient. 

The  paper  is  organised  as  follows.  In  §2  the  method  is  described  and 
the  gradient  estimator  is  constructed.  In  §3  the  convergence  properties  are 
considered.  Finally  in  §4  we  discuss  some  questions  related  to  the  numerical 
aspects  of  the  problem  and  we  provide  some  numerical  examples. 

§2.  The  Method 

Suppose  we  are  given  a set  of  scattered  points  in  a domain  Q C IR2 
S = {P*(xi,yi)\P*  €Q,i  = l,...,N}, 
and  the  set  of  functional  data 

F = {(P?,f(Pn),i  = 1 N}, 

where  f(x,  y)  is  an  unknown  function  defined  on  Q.  Without  loss  of  generality, 
we  suppose  Q = [0, 1]  x [0,1], 

The  first  step  of  our  method  consists  in  generating,  from  F,  a new  set  of 
functional  values,  say  F,  located  on  a dyadic  lattice  T of  Q, 

T = {Pk,  k = (/c1,fc2)€K2,  k = l,...,2fi}, 


that  is 

f = rnJih)),  PkGT}. 

It  is  clear  that  the  new  values  are  affected  by  errors 

f(Pk)  = f(Pk)  + e(Pk). 

Therefore  we  need  to  construct  F by  an  efficient  computational  method  such 
that  the  error  is  less  or  of  the  same  order  as  that  generated  by  the  derivative 
approximation  we  will  use.  A possible  strategy  is  to  interpolate  the  data  by  a 
local  method  of  a suitable  order  m,m  < a (for  instance  a moving  least  squares 
technique,  [6,  8])  which  gives  a smooth  interpolating  function  f(x , y)  £ Cm(Q) 
with  an  error  e(x,  y)  = 0(hm),  where  h is  a local  parameter  depending  on  the 
distribution  of  the  points  P?  in  Q (usually  h = 1 /VN). 

The  new  set  F can  be  thought  of  as  a sample  coming  from  a stochastic 
process  depending  on  the  points  P‘  £ S.  Then  we  can  use  the  method  de- 
scribed in  [4]  for  noisy  data.  In  the  next  section,  we  will  perform  a wavelet 
decomposition  in  order  to  smooth  the  errors  e,  and  we  will  define  an  estimator 
9j(n)(x,  y)  of  the  underlying  function  f(x,y).  Then  we  will  approximate  the 
unknown  gradient  using  the  function  estimator  gj^(x,y),  and  the  classical 
approach  of  central  finite  differences. 
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In  this  paper  we  assume  that 

i)  f{xi  V ) belongs  to  an  Holder  space  of  order  a > 2 on  Q*  D Q,  say  Ca(Q*  ). 

ii)  We  consider  a multiresolution  analysis  of  L2(1R2)  given  by  the  tensor 
product  of  two  one-dimensional  s-regular  multiresolution  analysis  such 
that  the  scaling  function  cj)  is  a coiflet,  that  is  a compactly  supported 
orthonormal  function  such  that  the  scaling  function  <j>( x)  and  the  wavelet 
ip(x)  have  L — 1 and  L,  vanishing  moments  respectively.  Moreover,  we 
assume  that  L > [a]  + 1 and  s > a. 

2.1.  The  function  estimator 

As  mentioned  above,  for  approximating  the  gradient  of  f(x,y),  we  will  use 
the  function  gj^){x>y)  defined  as  in  [4],  In  this  section  we  briefly  describe  its 
definition  and  the  motivations  that  lead  to  consider  this  function. 

It  is  known  that  when  we  perform  a MRA  using  data  given  on  a subset  of 
IR2,  we  need  to  take  into  account  the  problem  of  reducing  the  boundary  errors. 
To  this  end,  we  consider  a function  g(x,y)  £ C“(1R2)  compactly  supported 
on  Q*  D Q such  that 


g(x,y)  = fix,y),  V(x,r/)e<2, 

and  a new  dyadic  lattice  T*  on  Q*  of  dimension  2"  x 2",  such  that  T*  — 
{T  U { points  sampled  in  Q*  \ Q}}.  Moreover,  the  advantage  of  the  nested 
structure  of  a MRA  is  that  to  provide  an  efficient  tree-structure  algorithm 
for  the  decomposition  of  functions  in  V„  for  which  the  smoothing  coefficients 
are  given. 

In  applications,  a function  is  given  in  sampled  form,  and  it  is  therefore 
necessary  to  approximate  the  projection  Pyn  on  the  space  Vn,  by  some  oper- 
ator nn,  and  to  derive  a reasonable  estimator  of  nn  in  terms  of  the  sampled 
values.  The  choice  of  Hn  and  of  its  estimator  is  suggested  by  the  following 
facts  (see  [1,4]). 

The  set  of  nonzero  coefficients  (g,$n,k)  has  cardinality  equivalent  to 
0(22n).  Moreover, 

\(g,  $„,*>  - 2~ng(Pk)\  < C2~n2~na,  (1) 

where  C is  a constant  depending  on  the  smoothing  function  $(x,  y)  and  on 

g{x,y)- 

As  a consequence,  we  define 

(n„ff)(a:,  y)  = 2-n  ^ g(Pk)$n,k{x,y).  (2) 

pfc6r- 

Since  we  have  data  corrupted  by  the  interpolation  errors,  we  consider  the 
estimator  of  n„ 

(flng)(x,y)  = 2~n{  f(Pk)$n,k(x’V)  + iKpk)$n,k(a:,  y)}-  (3) 

Pk€T  Pk€T*\T 
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This  estimator  may  lead  to  an  oscillatory  solution  bearing  too  much  fidelity  to 
the  data.  Since,  for  numerical  differentiation,  we  need  to  correctly  smooth  the 
data,  we  have  to  associate  to  each  sample  of  size  2"  x 2n  a resolution  j(n)  < n, 
and  to  consider  the  orthogonal  projection  of  (nnff)(a:,  y)  onto  that  is 

&(»)(*.  3/)  = {PvjW^n9){x,y).  (4) 

The  parameter  j(n)  governs  the  smoothness  of  our  estimator,  and  it  is 
important  to  choose  it  in  the  right  way  because  it  controls  the  tradeoff  between 
the  fidelity  to  the  data  and  the  smoothness  of  the  resulting  solution.  From  a 
theoretical  point  of  view,  the  smoothing  parameter  must  tend  to  infinity  at 
the  correct  rate,  as  the  amount  of  information  in  the  data  grows  to  infinity. 

2.2.  The  gradient  estimator 

We  now  consider  the  construction  of  the  gradient  estimator.  When  dealing 
with  gridded  data,  it  is  natural  to  approximate  the  gradient  using  the  usual 
centered  difference  formulas.  Let 

(grad  f)(x,y)  = (/Cl)(x,  y),  f{v)(x,  2/))  , 

be  the  gradient  of  f(x,  y),  and  let  D*,  D'j  be  the  centered  difference  operators 
which  use  r equispaced  points  in  the  x or  y direction  respectively.  If  r < [a] , 
we  know  that 


(D\f)(Pk)  = f{t)(Pk)  + 0{ 2-"<r-1>),  t = [x,  y}. 


Then,  at  each  point  of  the  lattice  T,  we  approximate  (grad  f)(Pk)  by 


(grad  f)(Py)  = ((D*rgj{n))(Pk),  (D?ffm)(A))  ■ (5) 

Using  the  data  (5),  it  is  possible  to  define,  at  each  point  of  Q,  an  estimator  of 
f{t)(x,y),  t = {x,y},  by  the  operator  nn.  Namely,  we  consider  the  function 
g(x,y),  and  we  define 

(grad  f)(x,y)  = ((iW^Xa:,  y),  (flng^^x,  y)j  , (6) 


(iln9U)(x,y)  = 2~n{  £ (0&(*))(ft)*„,k(*,ff) 
Pk€T 

+ E S(t)W*  n,k(x,y)}. 
Pk€T"\T 


where 
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§3.  Asymptotic  Properties 

Under  the  assumptions  of  Section  2,  and  taking  into  account  the  properties 
of  the  projection  Pv, , we  know  that 

\(Pv,g)(x,y)  - g{x,y)\  = 0(2~la),  Vfr.t/jgQ*.  (7) 

Moreover,  from  (1)  we  have 

l(IItf)(*,lO  - (JV,s)(*,»)|  = 0(2-'“)  V(x,»)  € Q\  (8) 

then 

|(njff)(x, y)  - f(x,  y)|  = 0(2“'“)  V(x,  y)  £ Q.  (9) 

Using  relations  (7),  (8),  (9)  and  the  results  stated  in  [4],  we  have  proved 
Proposition  1.  If  assumptions  i)  and  ii)  of  Section  2 hold,  we  have 

\gm(x,y)  - f(x,y)\  = 0( 2"*fi>“)  + 0(hm), 

for  every  (x,  y)  £ Q. 

Remark  1.  This  result  points  out  how  the  choice  of  j(n)  depends  on  the 
sample  dimension  N and  on  m.  In  fact  it  has  to  be  chosen  so  that  2~^n'>ot  is 
less  or  of  the  same  order  of  hm. 

Proposition  2.  If  assumptions  i),  ii)  of  Section  2 hold,  the  approximation 
(6)  of  the  gradient  satisfies,  asymptotically,  the  following  bound 

|(ftnffw)(*,y)  - f(t)(x,y)  I =0(2-J<*>°)  + 0(2 -"(r-1)) 

+0(hm)  + 0(  2-"(“-1)), 


for  every  (x,  y)  £ Q. 


§4.  Numerical  Results 

In  this  section,  we  discuss  some  questions  related  to  the  computational  costs 
and  to  the  numerical  implementation  of  the  method  we  have  studied.  We  also 
present  some  numerical  results. 

4.1.  Computational  costs 

The  computational  costs  are  essentially  given  by  the  wavelet  decomposition 
and  by  the  construction  of  the  gridded  data  set  F of  dimension  22n.  For 
the  wavelet  decomposition,  they  are  at  the  most  of  the  same  order  of  the 
sample  dimension,  that  is  0(22n).  For  the  construction  of  F,  if  we  use  a local 
method  of  order  m N,  they  are  given  by  the  solution  of  N linear  systems 
of  dimension  2m  + 1.  Then  the  computational  costs  will  be  0((  2m3+1  )3N)  + 
0(22n). 
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4.2.  Numerical  implementation 

In  this  section  we  discuss  some  questions  related  to  the  construction  of 
9j(n){x,y)  and  (flng^)(x,y)  . Following  the  idea  of  §3.1,  we  need  to  extend 
the  given  signal  to  a suitable  square  Q*  D Q 

Q*  = [—2 ~nK,  1 + 2 ~nK]  x [-2 ~nK,  1 + 2 ~*K], 

forcing  it  to  be  zero  at  the  Q*  boundary.  This  is  necessary  in  order  to  avoid 
undesirable  behaviour  at  the  boundary.  The  extension  can  be  done  following 
a method  proposed  in  [2],  Note  that  K is  related  to  the  number  of  points  we 
consider  outside  Q.  On  one  hand,  it  cannot  be  chosen  too  small,  otherwise 
undesirable  boundary  oscillations  could  occur.  On  the  other  hand,  it  depends 
on  the  resolution  level  jf(n).  In  fact,  it  has  to  be  chosen  so  that  the  discrete 
wavelet  decomposition  can  be  performed. 

In  the  numerical  examples  we  have  used  the  coiflets  with  L = 6 vanishing 
moments.  For  the  pointwise  gradient  approximation,  we  have  used  the  central 
finite  difference  of  order  4 (r  = 5). 


4.3.  Numerical  examples 

In  this  paper  we  present  the  results  achieved  for  the  test  functions 


fi(x,y)  = 0.75 exp[—((9a;  - 2)2  + (9 y - 2)2)/4] 

+ 0.75  exp[— ((9x  + l)2/49  + (9  y + 1)2/10)] 
+ 0.5exp[-((9®  - 7)2  + (9 y - 3)2)/4] 

— 0.2exp[— ((9x  — 4)2  + {9y  — 7)2)], 

h{x,y)  = i 1 - ==, 

0(1  + 2exp(— 2^/lOCte2  + 100 y1  - 6.7) 


defined  on  the  unit  square  [0, 1]  x [0, 1].  We  have  considered  N scattered  points 
on  Q,  and  have  constructed  a new  gridded  data  set  F of  size  2"  x 2n  using 
the  modified  quadratic  Shepard  method.  The  smoothing  parameter  j(n)  has 
been  chosen  taking  into  account  Remark  1 of  Section  3.  Therefore,  having 
used  the  modified  quadratic  Shepard  method  with  n — 5,  a possible  choice  is 
j(fi)  = 4. 

We  now  present  our  results,  and  compare  them  with  those  obtained  with 
the  method  (L-method)  proposed  in  [8]  which,  among  those  we  find  in  the 
literature,  we  belive  is  preferable  both  for  its  theoretical  aspects  and  numerical 
performance. 

The  following  examples  show  how  the  proposed  method  provides  an  ap- 
proximation which  seems  to  have  the  same  behaviour  for  the  functions  con- 
sidered, both  for  graphical  results  and  for  errors.  For  the  sake  of  brevity,  we 
show  only  the  approximations  of  one  gradient  component,  but  give  the  relative 
errors  for  both  of  them. 
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Fig.  1.  Example  1,  TV  = 100.  On  the  left-hand  side,  the  result  of  our  method. 
On  the  right-hand  side,  the  result  of  the  L-method. 


Fig.  2.  Example  1.  N = 300.  On  the  left-hand  side,  the  results  of  our  method. 
On  the  right-hand  side,  the  result  of  the  L-method. 


Example  1.  Consider  N = 100  and  N — 300  values  of  fi(x, y).  In  Figs.  1 
and  2 we  present  the  ^-partial  derivative  approximations.  The  following  table 
lists  er—  relative  error  of  our  method  and  erL:=  relative  error  of  the  L-method: 


fl 

N = 

= 100 

TV  = 

= 300 

er 

erL 

er 

erL 

fx 

24% 

55% 

6% 

12% 

fv 

43% 

226% 

36% 

36% 

92 


M.  Bozzini  and  M.  Rossini 


Pig.  3.  Example  2.  N = 100.  On  the  left-hand  side,  the  result  of  our  method. 
On  the  right-hand  side,  the  result  of  the  L-method. 


Fig.  4.  Example  2.  N = 300.  On  the  left-hand  side,  the  result  of  our  method. 
On  the  right-hand  side,  the  result  of  the  L-method. 


Example  2.  Consider  N = 100  and  N = 300  values  of  y).  In  Figs.  3 
and  4 we  present  the  results  achieved  for  the  2/-partial  derivative,  and  in  Fig.  5 
the  error  functions  for  N = 300.  The  following  table  lists  er:=  relative  error 
of  our  method,  and  erL:=  relative  error  of  the  L-method: 


h N = 100  N = 300 

er  erL  er  erL 

35%  400%  6%  54% 

27%  64%  6%  135% 


r 

jv 


On  a Method  of  Numerical  Differentiation 


93 


Fig.  5.  Example  2.  N = 300.  On  the  left-hand  side,  the  error  function  of  our 
method.  On  the  right-hand  side,  the  error  function  of  L-method. 


0 o 


Fig.  6.  Example  2.  N = 300.  The  approximation  of  fc(x,  y).  The  absolute 
maximum  error  is  0.015. 

Finally,  as  is  usual  in  the  literature,  we  consider  an  application  to  function 
recovery  which  shows  the  goodness  of  the  gradient  approximation.  We  recover 
f2{x,y)  by  Hermite  interpolation  at  16  x 16  nodes,  where  we  interpolate  the 
data  coming  from  the  estimator  gj^(x,y)  and  from  the  approximated  gradi- 
ent (Fig.  6). 
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